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ABSTRACT 


The  interpretation  of  satellite-  or  aircraft-borne  synthetic  aperture  radar  (SAR) 
imagery  requires  increased  understanding  of  processes  that  change  the  state  of  the  sea 
surface.  Wind  stress  that  results  from  kilometer-scale  boundary  layer  circulations 
produces  sea-surface  stress  patterns  responsible  for  some  of  this  sea  surface  variability.  A 
Galerkin  model  of  these  boundary  layer  circulations — that  take  the  form  of  two- 
dimensional  rolls  and  three-dimensional  convective  cells — is  developed  here  to  study  their 
effect  on  sea  surface  stress. 

The  objectives  of  previous  investigations  using  spectral  models  were  focused  on  the 
structure  of  the  circulations  in  the  middle  levels  of  the  boundary  layer,  and  so  for 
simplicity  the  stress  and  heat  flux  were  assumed  to  vanish  at  the  lower  boundary. 

However  these  boundary  conditions  do  not  allow  the  study  of  stress  variability  at  the  sea 
surface.  The  goal  of  this  thesis  is  to  create  a  marine  atmospheric  boundary  layer  model 
that  allows  nonzero  heat  and  momentum  fluxes  at  the  lower  boundary.  Surface  layer 
similarity  theory  is  used  to  specify  constant  forcing  parameter  values  in  the  lower 
boundary  conditions.  The  focus  of  this  study  is  to  determine  whether  these  lower 
boundary  conditions  allow  the  model  to  capture  the  behavior  of  the  roll  circulations  and,  if 
they  do,  to  determine  the  fluxes  at  the  lower  boundary  that  influence  kilometer-scale  sea 
stress  variability. 

Case  study  results  demonstrate  that  the  model  does  indeed  capture  the  spatial 
organization  of  the  roll  circulations  and  does  correctly  reproduce  the  vertical  profiles  of 
heat  and  momentum  flux.  The  results  also  show  that  the  temperature  perturbation  is  a 
maximum  at  the  lower  boundary  where  the  primary  heat  source  is  located.  The  temporal 
periodicity  of  the  system  is  consistent  with  physical  and  numerical  analyses  as  well. 


However  there  are  two  distinct  case  study  results  that  indicate  further  redesign  of 
the  model  is  needed  before  the  surface  stress  patterns  can  be  studied.  The  first  result  is 
that  only  an  amplifying  roll  solution  could  be  found,  indicating  that  a  source  of  energy  is 
not  being  damped  properly.  The  other  result  is  that  the  maximum  vertical  velocity 
magnitude  occurs  at  the  lower  boundary,  which  produces  a  profile  that  is  not  consistent 
with  observations. 

The  case  study  results  point  to  the  need  to  re-evaluate  the  use  of  surface  layer 
similarity  theory  for  obtaining  the  parameter  values  in  the  new  boundary  conditions.  For 
example,  because  the  roll  circulations  scale  with  the  boundary  layer  depth,  it  seems  likely 
that  convective  scaling  theory  might  produce  more  appropriate  lower  boundary  conditions 
than  would  surface  layer  theory. 
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The  interpretation  of  satellite-  or  aircraft-borne  synthetic  aperture  radar  (SAR) 
imagery  requires  increased  understanding  of  processes  that  change  the  state  of  the  sea 
surface.  Wind  stress  that  results  from  kilometer-scale  boundary  layer  circulations 
produces  sea-surface  stress  patterns  responsible  for  some  of  this  sea  surface  variability.  A 
Galerkin  model  of  these  boundary  layer  circulations — that  take  the  form  of  two- 
dimensional  rolls  and  three-dimensional  convective  cells — is  developed  here  to  study  their 
effect  on  sea  surface  stress. 

The  objectives  of  previous  investigations  using  spectral  models  were  focused  on  the 
structure  of  the  circulations  in  the  middle  levels  of  the  boundary  layer,  and  so  for 
simplicity  the  stress  and  heat  flux  were  assumed  to  vanish  at  the  lower  boundary. 

However  these  boundary  conditions  do  not  allow  the  study  of  stress  variability  at  the  sea 
surface.  The  goal  of  this  thesis  is  to  create  a  marine  atmospheric  boundary  layer  model 
that  allows  nonzero  heat  and  momentum  fluxes  at  the  lower  boundary.  Surface  layer 
similarity  theory  is  used  to  specify  constant  forcing  parameter  values  in  the  lower 
boundary  conditions.  The  focus  of  this  study  is  to  determine  whether  these  lower 
boundary  conditions  allow  the  model  to  capture  the  behavior  of  the  roll  circulations  and,  if 
they  do,  to  determine  the  fluxes  at  the  lower  boundary  that  influence  kilometer-scale  sea 
stress  variability. 

Case  study  results  demonstrate  that  the  model  does  indeed  capture  the  spatial 
organization  of  the  roll  circulations  and  does  correctly  reproduce  the  vertical  profiles  of 
heat  and  momentum  flux.  The  results  also  show  that  the  temperature  perturbation  is  a 
maximum  at  the  lower  boundary  where  the  primary  heat  source  is  located.  The  temporal 
periodicity  of  the  system  is  consistent  with  physical  and  numerical  analyses  as  well. 
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However  there  are  two  distinct  case  study  results  that  indicate  further  redesign  of 
the  model  is  needed  before  the  surface  stress  patterns  can  be  studied.  The  first  result  is 
that  only  an  amplifying  roll  solution  could  be  found,  indicating  that  a  source  of  energy  is 
not  being  damped  properly.  The  other  result  is  that  the  maximum  vertical  velocity 
magnitude  occurs  at  the  lower  boundary,  which  produces  a  profile  that  is  not  consistent 
with  observations. 

The  case  study  results  point  to  the  need  to  re-evaluate  the  use  of  surface  layer 
similarity  theory  for  obtaining  the  parameter  values  in  the  new  boundary  conditions.  For 
example,  because  the  roll  circulations  scale  with  the  boundary  layer  depth,  it  seems  likely 
that  convective  scaling  theory  might  produce  more  appropriate  lower  boundary  conditions 
than  would  surface  layer  theory. 
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The  interpretation  of  satellite-  or  aircraft-borne  synthetic  aperture  radar  (SAR) 
imagery  requires  increased  understanding  of  processes  that  change  the  state  of  the  sea 
surface.  Wind  stress  that  results  from  kilometer-scale  boundary  layer  circulations 
produces  sea-surface  stress  patterns  responsible  for  some  of  this  sea  surface  variability.  A 
Galerkin  model  of  these  boundary  layer  circulations — that  take  the  form  of  two- 
dimensional  rolls  and  three-dimensional  convective  cells — is  developed  here  to  study  their 
effect  on  sea  surface  stress. 

The  objectives  of  previous  investigations  using  spectral  models  were  focused  on  the 
structure  of  the  circulations  in  the  middle  levels  of  the  boundary  layer,  and  so  for 
simplicity  the  stress  and  heat  flux  were  assumed  to  vanish  at  the  lower  boundary. 

However  these  boundary  conditions  do  not  allow  the  study  of  stress  variability  at  the  sea 
surface.  The  goal  of  this  thesis  is  to  create  a  marine  atmospheric  boundary  layer  model 
that  allows  nonzero  heat  and  momentum  fluxes  at  the  lower  boundary.  Surface  layer 
similarity  theory  is  used  to  specify  constant  forcing  parameter  values  in  the  lower 
boundary  conditions.  The  focus  of  this  study  is  to  determine  whether  these  lower 
boundary  conditions  allow  the  model  to  capture  the  behavior  of  the  roll  circulations  and,  if 
they  do,  to  determine  the  fluxes  at  the  lower  boundary  that  influence  kilometer-scale  sea 
stress  variability. 

Case  study  results  demonstrate  that  the  model  does  indeed  capture  the  spatial 
organization  of  the  roll  circulations  and  does  correctly  reproduce  the  vertical  profiles  of 
heat  and  momentum  flux.  The  results  also  show  that  the  temperature  perturbation  is  a 
maximum  at  the  lower  boundary  where  the  primary  heat  source  is  located.  The  temporal 
periodicity  of  the  system  is  consistent  with  physical  and  numerical  analyses  as  well. 
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However  there  are  two  distinct  case  study  results  that  indicate  further  redesign  of 
the  model  is  needed  before  the  surface  stress  patterns  can  be  studied.  The  first  result  is 
that  only  an  amplifying  roll  solution  could  be  found,  indicating  that  a  source  of  energy  is 
not  being  damped  properly.  The  other  result  is  that  the  maximum  vertical  velocity 
magnitude  occurs  at  the  lower  boundary,  which  produces  a  profile  that  is  not  consistent 
with  observations. 

The  case  study  results  point  to  the  need  to  re-evaluate  the  use  of  surface  layer 
similarity  theory  for  obtaining  the  parameter  values  in  the  new  boundary  conditions.  For 
example,  because  the  roll  circulations  scale  with  the  boundary  layer  depth,  it  seems  likely 
that  convective  scaling  theory  might  produce  more  appropriate  lower  boundary  conditions 
than  would  surface  layer  theory. 
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The  use  of  satellite-  or  aircraft-borne  synthetic  aperture  radars  (SAR)  as 
meteorological  and  oceanographic  observation  systems  holds  great  promise  (Visecky  and 
Stewart  1982).  Unfortunately,  the  interpretation  of  SAR  signals  is  greatly  dependent  upon 
the  nature  of  the  constantly  changing  state  of  the  sea  surface.  The  dominant  atmospheric 
mechanism  for  forcing  the  sea  surface  is  wind  stress,  and  so  a  better  understanding  of  the 
effects  of  this  stress  is  required  for  us  to  be  able  to  interpret  correctly  the  SAR  signals. 

One  approach  is  to  use  atmospheric  modeling  of  boundary  layer  flow  to  analyze  how  it 
might  modulate  sea-surface  stress  patterns. 

We  focus  our  study  on  kilometer-scale  marine  atmospheric  boundary  layer 
circulations.  There  are  two  types  of  kilometer-scale  circulations  that  affect  sea-surface 
stress.  The  first  type  comprises  quasi-two-dimensional  boundary  layer  rolls.  Evidence 
from  SAR  imagery  of  stress  patterns  caused  by  roll  circulations  is  well  documented,  for 
example,  by  Gerling  (1985,  1986).  The  second  type  comprises  convective  cells,  which  are 
fully  three-dimensional  flows  that  may  form  independently  of  the  rolls  in  less  windy 
conditions.  Such  cells  may  also  evolve  from  two-dimensional  rolls  as  the  forcing  rates 
increase.  As  described  extensively  elsewhere  (e  g.  Brown  1980;  Stensrud  1987),  each  of 
these  boundary  layer  circulations  may  be  driven  by  both  convective  and  dynamic  instability 
mechanisms.  Some  important  mechanisms  include  boundary  layer  wind  shear,  air/sea 
temperature  difference,  and  the  magnitude  of  heating  effects  in  clouds.  All  of  these 
mechanisms  are  ultimately  responsible  for  the  downward  transport  of  horizontal 
momentum  that  produces  stress  variability  at  the  sea  surface. 
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A  number  of  nonlinear  modeling  studies  of  boundary  layer  rolls  (Shirer  1986; 
Stensrud  and  Shirer  1988;  Laufersweiler  and  Shirer  1989;  Haack  and  Shirer  1992)  have 
sought  to  identify  the  dominant  forcing  mechanisms  and  characteristic  spatial  and  temporal 
responses  of  these  convective  boundary  layer  circulations.  These  nonlinear,  Galerkin 
models,  which  are  based  on  the  Boussinesq  system  of  equations  and  use  sinusoidal  basis 
functions,  provide  wind  profile  modifications,  preferred  orientation  angles  and 
wavelengths,  and  simplicity  in  the  study  of  bifurcation  and  stability  properties.  These 
models  also  yield  profiles  of  the  vertical  fluxes  of  horizontal  momentum  and  heat. 

The  objectives  of  the  above  studies  are  focused  on  determining  the  structure  of  the 
circulations  in  the  middle  levels  of  the  boundary  layer,  and  so  for  simplicity  the  stress  and 
heat  flux  are  assumed  to  vanish  at  the  lower  boundary.  However,  these  boundary 
conditions  cannot  allow  us  to  study  the  stress  variability  at  the  sea  surface,  unless  bulk 
aerodynamic  techniques  are  used  (Stull  1988).  Thus  the  goal  of  this  thesis  is  to  create  a 
nonlinear,  Boussinesq,  three-dimensional  boundary  layer  model  that  allows  nonzero  heat 
and  momentum  fluxes  at  the  lower  boundary.  There  are  a  number  of  approaches  that 
might  be  adopted  to  determine  appropriate  lower  boundary  conditions.  For  example, 
Schmidt  and  Schumann  (1989)  used  surface  layer  similarity  theory  to  specify  the  lower 
boundary  conditions  in  their  model,  which  they  used  to  successfully  capture  kilometer- 
scale  boundary  layer  structures. 

We  also  use  similarity  theory,  but  adopt  a  different  approach  to  specify  our 
boundary  conditions.  We  set  the  bottom  of  our  domain  at  10  m  above  the  sea  surface  and 
then  choose  lower  boundary  conditions  that  represent  wind  and  temperature  values 
typically  observed  at  this  height  in  the  surface  layer.  To  test  this  approach,  we  use 
empirically  derived  similarity  relationships  for  an  unstable  environment  (Stull  1988). 

These  relationships  lead  to  constant  forcing  parameters  that  may  be  used  to  specify  the 
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lower  boundary  conditions.  These  boundary  conditions  allow  for  a  more  realistic 
incorporation  of  thermal  forcing,  which  in  earlier  models  is  spread  throughout  the  entire 
vertical  domain  via  the  Rayleigh  number  Ra.  In  reality,  the  thermal  forcing  at  the  lower 
boundary  should  drive  the  system  and  the  value  of  the  traditional  Rayleigh  number  should 
be  inferred  from  the  model  solutions.  Finally,  and  for  simplicity,  the  upper  boundary 
conditions  are  taken  to  be  rigid,  stress-free,  and  perfectly  conducting  as  in  the  earlier 
models. 

We  then  develop  the  vertical  basis  functions  for  the  model  expansions  by  solving  a 
more  realistic  eigenfunction  problem  based  on  the  new  boundary  conditions.  These  basis 
functions  are  expressed  in  terms  of  exponential  and  trigonometric  functions.  When  we 
create  the  spectral  model,  products  of  these  basis  functions  must  be  integrated;  we  obtain 
expressions  for  these  integrals  using  the  symbolic  manipulator  software  DERIVE.  These 
results  are  then  incorporated  into  the  FORTRAN  program  that  does  the  numerical 
integration  of  the  model  and  subsequent  simulation  of  the  boundary  layer  circulations. 

The  focus  of  this  study  is  to  determine  whether  the  lower  boundary  conditions 
specified  using  similarity  theory  do,  in  fact,  allow  the  model  to  capture  the  behavior  of  the 
roll  circulations  and,  if  they  do,  to  determine  the  fluxes  at  the  lower  boundary  that 
influence  kilometer-scale  sea  stress  variability.  In  Chapter  2,  we  explain  the  development 
of  the  60-coefficient  nonlinear  dynamic  model,  which  is  still  much  smaller  than  a  Large 
Eddy  Simulation  (LES)  model,  and  so,  in  principle,  can  be  studied  more  completely.  We 
first  determine  the  appropriate  partial  differential  system  of  equations,  next  form  the 
spectral  system,  and  finally  describe  the  computer  model.  In  Chapter  3,  we  present  the 
results  of  some  model  integrations,  first  specifying  the  appropriate  model  parameter  values 
and  initial  conditions  and  then  analyzing  the  data.  Finally,  in  Chapter  4,  we  evaluate  these 
model  results  and  suggest  areas  for  further  analysis  and  study. 


Chapter  2 

MODEL  DEVELOPMENT 
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The  convective  nature  of  boundary  layer  roll,  or  secondary,  circulations  allows  us  to 
choose  a  set  of  approximate  equations  in  which  certain  terms  from  the  complete  set  of 
equations  are  neglected.  These  approximate  equations  form  what  is  known  as  the  shallow 
Boussinesq  system.  The  roll  circulations  are  represented  by  perturbations  superimposed 
on  a  basic  state  that  is  time-independent,  isopycnic,  hydrostatic,  and  horizontally  moving 
(Shirer  1986).  We  assume  that  the  perturbations  possess  a  three-dimensional  structure, 

and  so  we  must  consider  both  roll-perpendicular  and  roll-parallel  variations.  The 
horizontal  domain  is  infinite  and  cyclically  continuous  at  x  =  0,.y  =  0  and  x  =  Lx,y  =  Ly, 

where  Lx  and  Ly  are  the  roll  wavelength  components.  The  vertical  domain  ranges  from 
z  =  hLB  to  z  =  zT ,  where  hw  is  the  height  of  the  lower  boundary  and  zT  is  the  cloud  top 
height  or  inversion  base.  Only  the  part  of  the  circulation  above  is  specifically 
modeled;  the  part  below  is  parameterized  with  surface  layer  similarity  theory.  In 
dimensionless  form  the  domain  is  defined  by  0  <  x*  <  1 ,  0  <  y*  <  1,  and  0  <  z*  <  1  A 
schematic  representation  of  the  domain  is  given  in  Figs.  2. 1-2.2. 

The  basic  state,  known  as  the  conductive  state  and  denoted  by  subscript  o,  is  given 
as 


v.  =  l/(r)i+F(z)j  (2.1) 

P.=P.  (22) 

T.(z)=T„-\T(z/zt)  (2.3) 

p,(z)  =  P,  -  P-gz,  (2.4) 


5 


AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

SEA  SURFACE 


Fig.  2.1.  A  schematic  cross-section  of  the  model  domain  in  the  x-z  plane. 
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Fig.  2.2.  A  schematic  cross-section  of  the  model  domain  in  the  y-z  plane. 
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where  subscript  oo  denotes  the  value  at  the  lower  boundary,  U(z)  and  V(z)  are  the 
components  of  the  mean  horizontal  wind,  and  &2T  is  a  constant  lapse  rate  defined  as  the 
temperature  difference  between  the  lower  and  upper  boundaries.  The  three-dimensional 
secondary  circulation,  denoted  by  primes,  is  given  by 


u'(x,y,z,t)=  u(x,y,zj)-U(z) 

(2.5) 

v'(x,y,z,t)=v(x,y,z,t)-V(z) 

(26) 

w'(x,y,z,t )  =  w(x,y,z,t) 

(2.7) 

p'(x,y,z,t)  =  p(x,y,z,t)-poo 

(28) 

T'(x,y,z,t)=  T(x,y,z,t)-T0(z) 

(2.9) 

p'(x,y,z,t)  =  p(x,y,z,t)  -  p,(z). 

(2.10) 

In  the  Boussinesq  system,  compression  is  considered  negligible.  Other  approximations 
include  writing  the  equation  of  state  in  a  normalized  form  and  linearizing  the  pressure 
gradient  force  (Shirer  1987). 

2. 1 .  The  Partial  Differential  Equations 

The  Boussinesq  system  of  equations  for  the  perturbation  temperature  T  and  the 
perturbation  velocity  vector  v'  takes  the  form  (e  g.  Shirer  1986) 

— +v»vr+v'«vr  =  ^H',+xv2r-oV4r  (2.1 1) 

dt  zT 

^L-+  v«Vv'  +  H''-^-  +  v'«Vv'  =  — —  Vp'-  Jkx  \'  +  g— k  +  vV2v'-£V4v'  (2.12) 

dt  dz  pm  L 
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V«v'  =  0,  (2.13) 

where  k  is  the  constant  eddy  thermometric  conductivity  and  v  is  the  constant  eddy 
viscosity.  We  have  added  additional  terms  to  (2. 1 l)-(2. 12)  with  dissipation  parameters  a 
and  e  to  control  the  growth  of  any  numerical  instabilities.  A  dimensionless  form  of  the 
equations  is  used  based  on  the  definitions 


»  =  («.  +  vj)  . 

ZT  "LB 

(2.14) 

.  w'k 

w  = - 

zr  ~^lb 

(2.15) 

x  =  x*Lx 

(2.16) 

y=y'Ly 

(2.17) 

z  —  z  (zT  — 

(2.18) 

t  _  f  iZT  ~hui) 

K 

(219) 

T,_  T  vkT„ 

S{ZT  ~  h LB  ) 

(2.20) 

p>=  pp~k1 

\ZT  ~  h LB  )4 

(2.21) 

C7(r)  =  |v(zr)|t/-(r-) 

(2.22) 

F(Z)  =  |V(i7^-(z-) 

(2.23) 

rs 

.V  J 

^  1 
«sj~^ 
II 

(2.24) 

g 


These  dimensionless  equations  take  the  form 


er 


d\“ 

dt 


-  =  V2r  -  y0V4r  -  Re  V*  •  vr*  +  Raw'  -  v*  •  VT 


(2.25) 


7-  =  -axVp*  +  PV2v*  -  //vv*  -f*Pk  x  v*  +  PFlk 


-(Re  V*  +  v*)«  Vv*  -  Re 

oz 


(2.26) 


where 


V=ax^ri+a^j+-^k.  (2.27) 

dx  dy  dz 

The  parameter  ft  is  the  dimensionless  version  of  a  and  the  parameter  p  is  the 
dimensionless  version  of  s .  The  values  of  the  dimensionless  parameters  fi  and  p  control 
the  magnitude  of  the  numerical  dissipation  terms.  The  dimensionless  forms  lead  to  two 
forcing  parameters  in  (2.25)-(2.26).  The  Reynolds  number  Re  is  given  by 

Re  =  |V(2r)|(zr -/,„)/*■  (2  28) 

and  represents  dynamic  forcing  by  the  mean  wind  shear  of  the  basic  state.  The  Rayleigh 
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number  Ra  is  given  by 


Ra  = 


gA,T(zT-hJ 

vkT^ 


(2.29) 


and  represents  thermodynamic  forcing  across  the  domain. 

Three  other  dimensionless  variables  appear  in  the  system.  The  eddy  Prandtl 
number  is  defined  as 


P  =  vJk. 

(2.30) 

The  roll  aspect  ratios,  ax  and  ay,  are  defined  as 

II 

i 

sf" 

(231) 

_  (*r  ^lb) 

y~  L 

y 

(2.32) 

The  vector  V*(/)  represents  the  dimensionless  mean  wind  profile. 

As  a  result  of  the 

definitions  (2.22)-(2.23),  it  follows  immediately  that  |V*(l)|  =  1. 


We  further  transform  the  momentum  equation  (2.26)  by  taking  the  curl  to  form  a 
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r 

f 

! 


i 

: 


vorticity  equation,  thereby  eliminating  pressure.  The  resulting  equation  is 


x  v*) 

= x  v*)  ~  x  v’)  “  x  (k  x  v*)] + x  rk) 


-Re[Vx(v**Vv*)]-Re 


Vx  w 


.d\'' 

dz* 


-[Vx(v*«Vv*)].  (2.33) 


In  order  to  solve  the  system  of  equations  (2.25)  and  (2.33),  we  must  specify 
appropriate  boundary  conditions.  For  simplicity,  the  lateral  boundaries  are  cyclic,  and  the 
upper  vertical  boundary  is  rigid,  stress-free,  and  perfectly  conducting.  In  contrast,  an 
important  source  of  energy  for  the  circulations  exists  at  the  lower  boundary;  there 
boundary  conditions  are  specified  using  surface  layer  similarity  theory  as  described  in 
Appendix  A.  The  vertical  boundary  conditions  are  summarized  below  for  Sj  >  0  and 

s»>0: 


r(°)+Sr^=° 

(2.34) 

•Wga- 

(2.35) 

v<°)-s>  J*)=0 

(2.36) 

dw\Q)  _  <?V( 0) 

dz  m  dz1 

(2.37) 

r(i)=o 

(2.38) 

(2.38) 
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du{  l)  <5V*(l) 
dz  ~  dz 


h-‘(1)  =  0. 


(2.39) 


where  sm,  the  Schramm  momentum  constant,  and  sT,  the  Schramm  temperature  constant, 
which  are  both  defined  in  Appendix  A,  control  the  forcing  at  the  lower  vertical  boundary. 

2.2.  The  Spectral  System 

The  next  step  in  the  development  of  the  model  is  to  create  an  appropriate  system  of 
ordinary  differential  equations  from  the  partial  differential  system  given  in  the  previous 
section.  We  use  the  Galerkin  technique  or  spectral  method  to  form  this  spectral  model. 
First,  we  substitute  the  Fourier  expansions  for  the  dependent  variables  into  the  appropriate 
partial  differential  equation.  Next,  we  multiply  the  resulting  equation  by  the  particular 
eigenfunction  that  the  amplitude  coefficient  multiplies  in  the  expansion.  We  then  integrate 
the  result  over  the  spatial  domain,  using,  where  possible,  orthogonality  of  the 
eigenfunctions  in  the  expansion.  Finally,  we  retain  only  those  terms  that  are  nonzero 
(Higgins  1987). 

Owing  to  the  complexity  of  the  model,  the  above  steps  are  performed  separately  for 
each  term  in  the  partial  differential  equations,  with  the  results  being  stored  in  an 
appropriate  matrix.  The  truncation  chosen,  which  is  listed  below,  results  in  a  system  of  60 
equations,  20  in  temperature  and  40  in  velocity.  Each  term  in  each  equation  may  be 
represented  by  a  60  x  60  matrix  that  is  multiplied  by  a  60-element  vector  representing  the 
60  amplitude  coefficients.  The  resulting  system  takes  the  form 


Ay  =  By  +  C(y)y, 


(2.40) 


where  A  is  a  60  x  60  matrix  of  inner  products  multiplying  the  60-element  vector  y  that 
contains  the  temporal  derivatives  of  the  amplitude  coefficients,  B  is  a  60  x  60  constant 
matrix,  based  on  the  linear  terms  of  the  partial  differential  equations,  multiplying  the  60- 
element  vector  y  that  contains  the  amplitude  coefficients,  and  C  is  a  60  x  60  ^-dependent 
linear  matrix,  based  on  the  advective  terms  of  the  partial  differential  equations,  multiplying 
the  same  60-element  vector.  Because  A  is  invertible,  we  are  then  able  to  solve  for  the 
time-dependent  amplitude  coefficients  by  numerically  integrating  the  following  system: 

y  =  A'By  +  A'C(y)y.  (2.41) 

To  find  (2.40)  from  (2.25)  and  (2.33),  we  first  choose  a  truncated  Fourier  expansion 
for  the  dependent  variable  T  composed  of  temporally  dependent  amplitudes,  spatially 
dependent  trigonometric  functions  with  horizontal  wavenumbers  /  and  m,  and  vertical 
basis  functions: 


r(x-y,z/) = it  r,  (<>«(*•, (2.42) 

1=0  j= 1 

where 

<"&(**>/)  =  1  (243) 

trigx (x*,y*)  =  sin  2rdx*  sin  limy  (2.44) 

trig2(x\y)  =  sin 2 jdx  cos2 nmy  (2.45) 

trig3(x\y*)  =  cos2;rfx*  sin  2  nmy*  (2.46) 
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*»&(**./)  =  cos2  jdx*  cos2wny* 

(2.47) 

Fj(z*)  =  coshm.z* - —  sinh<y,z*  for  Sj  <  1 

W 

(2.48) 

Fi(z*)  =  cosm,z* - —  sin<u,z*  for  sr>l 

(2.49) 

sTo)i  =  tanh/y,  for  s,.  <  1 

(2.50) 

sTd)]  =  tan o)x  for  Sj  >  1 

(2.51) 

F(z*)  =  cos<y,z* - —  sinrn  z*  for  j  =  2,3,4 

Sj-OJj 

(2.52) 

SyCOj  =  tan<a;  for  j  -  2,3,4, 

(2.53) 

and  the  vertical  basis  functions  are  orthogonal.  A  complete  discussion  of  the  vertical  basis 
functions  and  associated  eigenvalues  is  found  in  Appendix  B. 

We  next  must  choose  an  expansion  for  the  velocity  vector  v* .  Because  it  is 
nondivergent,  it  will  be  convenient  to  represent  v*  (nonuniquely)  as  the  curl  of  a  vector 
streamfunction  A*  that  is  the  sum  of  two  vector  streamfunctions  y/'  and  17*.  Therefore, 

the  velocity  perturbation  vector  v*  may  be  expressed  as 

v*  =  V  x  A*  =  V  x  ^*  +  V  x  17*.  (2.54) 

We  then  substitute  this  new  definition  for  v*  into  the  partial  differential  system  (2.25)  and 
(2.33).  We  can  now  use  a  truncated  Fourier  expansion  to  represent  the  dependent  vector 
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streamfunctions  yt*  and  rf : 

v'{x\y\z'/)=  XZ^«(/*Wp(Jf*>>',K(z*)i  +  0i  +  0k  (2  55> 

P=  0  «= 1 

i/*(x\/,z*/)=  0i  +  £ E^«('*)^(xVh(**)i  +  0k’  (2  56) 

P=Q  9=1 

where  trigr{x* ,y*)  is  defined  as  in  (2.43)-(2.47),  and 

hq{z )  =  cos  nr  ,2*  -  smmq  sin  mqz'  for  q  -  1,2,3,4  (2.57) 

smtaq  =  cot  or,  for  q  =  1,2,3,4,  (2.58) 

where  the  vertical  basis  functions  are  nonorthogonal.  Once  again,  a  complete  discussion 
of  the  vertical  basis  functions  and  eigenvalues  is  found  in  Appendix  B. 

Next,  we  substitute  the  temperature  expansion  (2.42)  and,  when  required,  the 
streamfunction  expansions  (2.55)-(2.56)  into  each  term  of  the  thermodynamic  equation 
(2.25).  We  then  multiply  each  term  by  the  appropriate  basis  function  and  integrate  over 
the  domain  0  x*  <  1, 0  <  >>*  <  1, 0  <  z*  <  1 .  We  use  the  following  identities  to  simplify  the 
equations: 


J  sin  lidx  sin  2 jdx'dx*  =  (2.59) 

0 

1 

J  sin  2  jdx*  cos2 ittx'dx  -  0 

0 


(2.60) 
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i 

J  sin  Ijony *  sin  Ijtmy  dy  =  y  (2.61) 

0 

1 

|  sin  2  jtmy  cos2 mny'dy *  =  0,  (2.62) 

0 


and  realize  that  only  integrals  involving  squares  of  horizontal  functions  are  nonzero.  We 
thus  conclude  for  the  linear  terms  that  we  only  obtain  nonzero  integrals  when  /  =  k  in  the 
equation 


|  J  rig,  (x\y)trigk  (x\y)±e*dy  =  y.  (2.63) 

0  0 

For  the  nonlinear  term  in  the  thermodynamic  equation,  we  must  consider  integrals  of  triple 
products  and,  therefore,  must  establish  selection  rules  to  determine  the  nonzero  integrals. 
These  rules  specify  that  the  integral 


jjrtgi(x\y*)pigk(x\y')trigm(x\y')dx'dy'  =  y,  (2.64) 

0  0 

exactly  when 

i  =  0  and  k  =  m 
or  k  =  0  and  i-m 

or  m-  0  and  i  =  k. 
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We  also  observe,  using  the  Lagrange  Identity  for  orthogonal  functions,  that 

=  y*/,  (2.65) 

0 

and  so  we  can  eliminate  expressions  containing  this  integral.  The  remaining  integrations 
of  vertical  basis  functions  are  evaluated  using  the  DERIVE  software  integration  routine 
and  simplifications  given  by  the  following  relationships. 

sT(0i  =  tanlwu,,  (2.66) 


sin©,  = 


cos<y;  = 


(i +5/(0  ;f2 

-(-i  y 

(1 


(267) 


(2.68) 


stq)j  =  tan  (Oj, 


(2.69) 


for  y  =  2,3,4, 


sin  mq 


-(-If 


and 


(2.70) 
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cos  mq 


(1  +  -s»2^2)/2 


smm9=cotmq. 


(2.71) 

(2.72) 


for  q  =  1, 2,3,4.  Integrations  of  vertical  basis  functions  in  the  mean  wind  advection  term 

i 

are  of  the  form  J  V*(z*]Fj(z')Fl(z*)dz*and  are  evaluated  using  an  IMSL  subroutine. 

0 

We  next  substitute  expansions  (2.55)-(2.56),  using  the  relationship  given  by  (2.54), 
into  each  term  of  the  vorticity  equation  (2.33).  We  once  again  multiply  each  term  by  the 
appropriate  basis  function  and  integrate  over  the  spatial  domain.  We  use  the 
trigonometric  integral  identities  (2.59)-(2.62)  to  eliminate  all  expressions  except  those  of 
the  form 


Wtrigpi*  S)tr'gXx' =  X’ 


0  0 


P  =  c- 


(2.73) 


For  the  nonlinear  term,  we  retain  expressions  of  the  form 


J  j ^iSP(x\y)trigXx\y)trigXx\y)dx,dy  =  y,  (2.74) 

0  0 


where 


p-  0  and  r  =  c 
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or  r  -  0  and  p-c 

or  c  =  0  and  p  =  r. 

Owing  to  the  nonorthogonality  of  the  momentum  vertical  basis  functions,  we  cannot 
eliminate  any  expressions  based  on  integrals  of  these  functions.  The  integrals  are 
evaluated  using  the  DERIVE  software  integration  routine  and  the  relationships  (2.66)- 

1 

(2.72).  The  mean  wind  advection  term  integrals  of  the  form  (2*^(2*)^’ are 

0 

also  evaluated  using  an  IMSL  subroutine;  integrals  involving  dV  jdz'  can  be  put  in  the 
above  form  via  integration  by  parts  (e  g.  Shirer  1986). 

2.3.  The  Computer  Model 

The  computer  code  used  to  represent  the  physical  and  mathematical  model 
presented  in  the  previous  sections  is  written  in  the  FORTRAN  programming  language.  It 

is  compiled  using  the  WATFOR77  compiler.  The  first  part  of  the  code  determines  the 
values  for  o)j  for  j  =  1,  2,  3, 4,  and  mq  for  q=  1,2,  3, 4,  based  on  input  values  for  the 

constants  sm  and  sT .  The  values  of  these  constants  are  provided  by  the  subroutine 
discussed  in  Appendix  A.  An  iterative  technique  is  used  to  solve  (2.50)-(2.51),  (2.53), 

and  (2.58)  for  the  eigenvalues.  Prior  to  integrating  the  system,  we  must  also  input  values 
for  the  parameters  /,  m,  az,ay,  p,  p,  Ra,  Re,  and  P. 

The  determined  and  input  values  are  used  to  fill  the  60  x  60  matrices  representing 
the  terms  in  the  spectral  system.  The  individual  matrices  that  compose  A,  B,  and  C  are 
grouped  as  temporal  derivative  term  matrices,  linear  term  matrices,  and  nonlinear  term 
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matrices.  We  can  rewrite  (2.25)  and  (2.33)  in  terms  of  this  matrix  notation: 

Ay  =  + Ct(y)y  (2.75) 

Aj=B1y+C1(y)y.  (2.76) 

All  the  matrices  in  a  given  group  are  then  added  together  . 

A=A]+A2  (2.77) 

B  =  B,+B2  (2.78) 

C(y)  =  Q{y)  +  C2{y).  (2.79) 

This  gives  the  system  (2.40)-(2.41).  The  temporal  derivative  term  matrix  and  the  linear 
term  matrix  are  filled  prior  to  the  integration  because  A  and  B  do  not  depend  on  the 
spectral  coefficient  in  y.  An  IMSL  routine  computes  A~x  and  a  different  IMSL  routine 

computes  A  'B.  The  nonlinear  term  matrix  must  be  filled  before  each  integration  time 
step  because  C  depends  on  y.  The  calculation  A  'C(y)  must  also  be  done  before  each 

time  step. 

The  IMSL  integration  routine  uses  the  Adams-Gear  method  and  requires  user- 
supplied  initial  conditions  for  the  60-element  vector  y.  At  each  time  step,  an  IMSL 
routine  calculates  A~]By  and  A~'C(y)y .  The  two  resulting  60-element  vectors  are  added 

together  to  give  the  right  sides  of  the  60  equations  for  which  we  are  solving.  The 
integration  routine  uses  these  results  to  calculate  new  values  for  the  60  time-dependent 
amplitude  coefficients  after  each  time  step.  We  integrate  the  model  until  the  square  root 
of  the  sum  of  the  squares  of  the  coefficients  equilibrates.  After  the  model  integration  is 
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complete,  these  amplitude  coefficients  are  used  to  construct  the  model  solution  using 
(2.4.  and  (2.55)-(2.56).  Examples  of  this  solution  are  given  in  the  following  chapter. 


Chapter  3 

MODEL  EVALUATION 
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In  order  to  investigate  whether  the  model  solutions  represent  the  expected  physical 
processes,  we  performed  a  series  of  model  integrations  and  analyses.  The  integrations 
included  varying  the  physical  parameter  values,  changing  the  initial  conditions,  or  using 
various  numerical  integration  techniques.  The  analyses  included  various  eigensystem 
evaluations  and  graphical  representations  of  the  evolution  of  the  solution.  These  various 
procedures  were  used  not  only  to  test  the  validity  of  the  physical  assumptions  made,  but 
also  to  identify  and  correct  any  computer  programming  errors.  The  results  of  these 
integrations  and  analyses  are  presented  in  Section  3.1. 

Once  we  chose  the  physical  parameter  values,  we  performed  a  simple  case  study  to 
evaluate  the  quality  of  the  model  solutions.  This  study  allows  us  to  examine  the  patterns 
produced  by  the  roll  circulations  and  the  associated  vertical  transfers  of  heat  and 
momentum.  By  comparing  these  results  with  those  of  other  models,  we  can  assess  the 
performance  of  our  model  and  the  effects  of  incorporating  the  forcing  rates  in  the  lower 
vertical  boundary  conditions.  The  results  of  this  case  study  are  given  in  Section  3.2. 

3.1.  Parameter  Values  and  Initial  Conditions 

The  first  physical  parameter  values  to  be  determined  are  those  given  by  the  basic 
states  of  the  atmosphere  and  sea  surface.  These  parameter  values  are  input  into  the 
subroutine  described  in  Appendix  A  that  determines  values  for  the  surface  roughness,  z0, 
and  the  Monin-Obhukov  length,  L.  The  values  obtained  for  z0  and  L  are  then  used  in 
another  subroutine  that  calculates  values  for  the  two  constant  lower  boundary  forcing 


parameters,  sm  and  sT .  The  parameter  values  used  in  the  model  case  study  are  given  in 
Table  3.1. 
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Table  3.1.  Parameter  values  used  in  the  model. 


PARAMETER 

VALUE 

ASPECT  RATIO  IN  THE  X-DIRECTION  (ax) 

0.001 

ASPECT  RATIO  IN  THE  Y-DIRECTION  ( ay ) 

0.7 

HEIGHT  OF  TOP  OF  DOMAIN  (ZT) 

1000  m 

HEIGHT  OF  BOTTOM  OF  DOMAIN  (hT  R) 

10  m 

U  MEAN  WIND  SPEED  AT  Z  =  ZT 

2.5  ms-1 

U  MEAN  WIND  SPEED  AT  Z  =  hT  R 

2.5  ms-1 

V  MEAN  WIND  SPEED  AT  Z  =  ZT 

2.5  ms'1 

V  MEAN  WIND  SPEED  AT  Z  =  hr  R 

2.5  ms*1 

GRAVITY  (x) 

9.8  ms*2 

CORIOLIS  PARAMETER  (/) 

0.0  s*1 

EDDY  VISCOSITY  (v) 

10  m2s*1 

EDDY  CONDUCTIVITY  ( k ) 

10  m2s*1 

PRANDTL  NUMBER  (P) 

1.0 

SEA  SURFACE  TEMPERATURE  (FJ 

293  K 

TEMPERATURE  AT  Z  =  hT  R 

5°  C 

MIXING  RATIO  AT  SEA  SURFACE 

6.0  g  kg*1 

MIXING  RATIO  AT  Z  =  hT  R 

15  g  kg*1 

RAYLEIGH  NUMBER  (Ra) 

2.75  x  10s 

REYNOLDS  NUMBER  (Re) 

350 
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In  an  effort  to  keep  the  system  simple  for  the  case  study,  we  use  a  constant  wind 
speed  of  2.5  ms "  in  both  the  x-  and  ^-directions.  This  gives  a  total  wind  vector  of 
magnitude  3.54  ms  ' .  We  choose  aspect  ratios  of  0.001  in  the  x-direction  and  0.7  in  the>>- 
direction  to  produce  a  quasi-two-dimensional  domain  oriented  approximately  in  the^-z 
plane.  The  most  important  parameter  values  supplied  to  the  initial  subroutine  are  those  for 
the  sea  surface  temperature  and  the  10  m  temperature.  The  surface  layer  vertical 
temperature  difference  AT  calculated  from  these  two  parameters  determines  the  amount 
of  thermal  energy  being  supplied  to  the  system  through  the  lower  boundary.  This  AT 
differs  from  AZT,  which  is  the  basic  state  temperature  difference  between  heights  zT  and 
hjj, .  In  Table  3 .2  we  present  various  values  of  AT  for  fixed  values  of  the  other 
parameters,  together  with  the  associated  values  of  z„,  L,  sm,  and  Sj. . 

Table  3.2.  Parameter  values  for  surface  layer  AT  and  associated  calculated  values  for  z, ,  L,  sm  and  sT . 


A  T(°C) 

L(m) 

Sm 

sT 

V 

0.27425x10" 

-4.6468 

0.277716 

0.491370 

5° 

0.27533  xlO-4 

0328111 

mam 

10° 

HM9 

■m 

BHIB 

1 

mm 

BH 

1.004940 

Bfi 

0.28134  x  10" 

BlPlifliB 

0.402118 

BIBB 

For  our  case  study,  we  choose  the  relatively  strong  forcing  value  of  AT  =  15°C  and  use 
the  associated  values  shown  in  Table  3 .2  for  sm  and  sT .  We  note  that  sT  >  1  here  so  that 
the  first  temperature  eigenfunction  is  given  by  (2.49)  and  (2.51).  Smaller  values  for  AT 
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would  require  us  to  use  the  hyperbolic  first  temperature  eigenfunction  given  by  (2.48)  and 
(2.50).  This  hyperbolic  case  remains  untested. 

Once  we  have  determined  the  values  for  sm  and  sT ,  the  computer  program,  using 
relationships  (2.51),  (2.53),  and  (2.58),  calculates  the  eigenvalues  <ay  and  mq  for 

j  =  q  =  1, 2, 3, 4.  These  values  for  the  case  study  are  given  in  Table  3.3.  We  can  use  these 
values  and  the  vertical  basis  functions  (2.49),  (2.52),  and  (2.57)  to  plot  profiles  for 
temperature  and  vertical  velocity.  These  plots  are  given  in  Figs.  3. 1-3.2.  They  indicate 
that  temperature  and  vertical  velocity  are  well  correlated.  The  temperature  profile  reveals 
that  the  eigenfunction  for  wavenumber  1  gives  a  linear  relationship,  as  seen  in 
observations,  that  is  consistent  with  the  basic  state,  and  succeeding  eigenfunctions 
superimpose  trigonometric  functions  on  this  linear  relationship.  The  maximum  in 
temperature  at  the  lower  boundary  is  consistent  with  the  thermal  forcing  occurring  there. 
The  vertical  velocity  profile  is  similar,  and  although  largest  at  z*  =  0,  a  maximum  in 
vertical  velocity  conceivably  could  occur  above  the  lower  boundary  through  interaction  of 
the  different  eigenfunctions. 


Table  3.3.  Case  study  values  for  <o)  and  m9  for  j  =  q  =  1, 2, 3, 4. 


(0  ■ 
j 

j=1  =  l 

0.121378 

1.153289 

j  =  q  =  2 

4,494503 

3.748059 

j  =  q  =  3 

7.725888 

6.655595 

j  =  q  =  4 

10.904572 

9.686992 

TEMPERATURE  BASIS  FUNCTION  (DIMENSIONLESS) 


Fig.  3.1.  Temperature  vertical  basis  function  profiles  for  the  four  wavenumbers  chosen. 


T 


T 


T 


-4.00  -2.00  0.00  2.00  4.00 

VELOCITY  BASIS  FUNCTION  (DIMENSIONLESS) 

Fig.  3.2.  Velocity  vertical  basis  function  (vertical  velocity)  profiles  for  the  four  wavenumbers  chosen. 
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The  initial  conditions  represented  by  the  initial  values  of  the  60  time-dependent 
amplitude  coefficients  are  determined  using  an  eigensystem  analysis.  We  find  the 
eigenvectors  associated  with  the  conjugate  pair  of  eigenvalues  whose  real  parts  change 
sign  with  a  change  in  the  forcing  rate.  Thus  these  eigenvectors  tell  us  the  relative 
magnitudes  of  the  coefficients  that  would  approximate  their  values  in  the  nonlinear 
branching  solution.  We  can  multiply  these  eigenvectors  by  a  constant  to  obtain  values  on 
the  scale  of  the  model  amplitude  coefficients  using  the  dimensionless  conversions,  and  so 
'  the  resulting  values  as  initial  conditions 

As  noted  in  Chapter  2,  the  IMSL  integration  routine  uses  the  Adams-Gear  method. 
A  series  of  model  integrations  were  performed  to  determine  the  proper  value  for  the  time 
step.  The  time  step  was  reduced  until  the  number  of  internal  subroutine  calls  was  small 
enough  to  allow  accurate  solutions  at  each  step.  We  determined  the  largest  possible  time 
step  to  be  0. 1  seconds,  which  corresponds  to  the  dimensionless  time  value  of  1.02  x  10-6. 
Model  integrations  were  run  for  5  hours,  which  corresponds  to  180,000  time  steps,  to 
allow  sufficient  time  for  the  system  to  equilibrate. 

Because  we  expect  the  roll  solutions  to  have  time  scales  on  the  order  of  an  hour  or 
less,  we  neglect  the  Coriolis  term  in  the  model.  We  specify  only  one  wavenumber  in  the  x- 
and  ^-directions,  and  so  we  set  /  =  1  and  m  =  1.  The  values  chosen  for  v  and  k 
determine  the  value  of  the  Prandtl  number  P  using  (2.30).  In  this  case,  v-k  =  10  mV , 
which  gives  P  =  1 .  We  also  choose  /?  =  n  =  I  for  the  numerical  dissipation  terms.  The 
value  for  the  dimensionless  Reynolds  number  Re  is  found  using  (2.28)  based  on  a  value  of 
the  wind  speed  (3.54  ms'1)  at  the  top  of  the  domain  and  a  domain  height  of  1000  m,  the 
result  is  Re  =  350.  The  value  used  for  the  dimensionless  parameter  Ra  is  determined 
based  on  a  stability  analysis. 
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Ideally,  we  would  choose  a  value  for  Ra  that  produces  a  conjugate  pair  of 
eigenvalues  having  slightly  positive  real  parts  in  the  linear  stability  analysis.  This  result 
signals  a  Hopf  bifurcation  producing  a  periodic  solution  that,  if  supercritical  and  hence 
stable,  is  realized  after  the  system  transients  have  disappeared  (e  g.  Pyle  1987). 
Unfortunately,  during  the  integrations  performed  for  this  study,  the  bifurcating  model 
so!  Ions  did  not  equilibrate,  indicating  a  source  of  energy  that  is  not  currently  being 
damped  properly.  This  problem  will  be  discussed  further  in  Chapter  4.  We  therefore 
choose  a  case  for  which  all  the  real  parts  of  the  eigenvalues  are  negative  and  so  the 
solution  is  damping.  By  finding  a  value  for  Ra  that  produces  a  very  slowly  damping 
solution,  we  can  simulate  a  system  approaching  an  attractor.  We  choose  Ra  =  2.75  x  105 
that  corresponds  to  A  J  =  0.85  °C  across  the  domain.  This  allows  us  to  integrate  the 
model  without  obtaining  excessive  growth  or  decay  on  time  scales  appropriate  for 
boundary  layer  rolls. 

By  plotting  the  square  root  of  the  sum  of  the  squares  of  the  60  time-dependent 
amplitude  coefficients,  we  can  determine  when  the  system  has  reached  a  quasi-equilibrated 
state.  A  plot  of  the  square  root  of  the  sum  of  the  squares  is  given  in  Fig.  3.3.  This  plot 
indicates  that  it  takes  approximately  5  hours  for  the  major  transients  to  disappear  using  the 
initial  conditions  determined  from  the  eigenvectors.  We  have  considered  the  possibility  of 
a  subcritically  branching  bifurcating  solution;  but  integrating  the  system  further  reveals 
that  the  solution  continues  to  damp  very  slowly,  indicating  that  such  a  subcritical  solution 
is  not  likely. 

We  can  also  checx  the  validity  of  the  model  by  examining  the  periodicity  of  the 
solution.  Using  a  linear  eigensystem  analysis,  we  can  determine  the  values  for  the 
imaginary  parts  Im  of  the  eigenvalues.  The  limiting  period  of  the  bifurcating  solution  is 


determined  by  calculating  (Pyle  1987) 


Period  =  —  (3  1) 

Im 

This  value  must  then  be  converted  to  dimensional  units  using  (2. 19).  Based  on  our 
analysis,  Im  =  1086.  Using  (3.1),  we  find  a  dimensionless  period  of  5.786  x  1 0" 3 ,  or  9.45 
minutes. 

Given  our  constant  wind  profile,  this  period  should  correspond  to  the  speed  of 
movement  of  the  rolls  caused  by  advection  by  the  mean  wind.  To  find  this  speed,  we  first 
note  that  the  ratio  Ly/ Lx  -  axjay  serves  as  a  measure  of  the  orientation  angle  0  of  a  roll- 
type  solution,  for  which  tan0  =  ax/ay ;  here  6  is  the  angle  between  the  roll  axis  and  east. 
At  any  level,  the  background  wind  angle  0  is  given  by  tan  $  =  V/U  =  V*/ U ’ .  Thus  it 
follows  that  the  angle  <f>-6  between  the  background  wind  vector  and  the  roll  axis  obeys 


tan(0  -  9) 


tan^-tanfl  _  ay*  ~axU* 
l  +  tan^tan0  ayU*+axV* 


(3.2) 


For  the  case  discussed  here,  ax  =  0.001  and  ay  =  0.7,  the  rolls  are  aligned  nearly  along  the 

x-axis,  and  so,  with  the  constant  values  U  =  V  =  2.5  ms' ,  the  angle  <f>-9  between  the 

mean  wind  direction  and  the  roll  is  approximately  45° .  This  result  is  also  given  by  (3.2)  in 
which  tan(^  -  0)  «  1  for  ax  «  ay  and  U*  =V\  This  means  that  advection  is  accomplished 

predominantly  by  the  V  component,  as  expected  from  (2.25)-(2.26)  when  az  «  ay. 

From  the  above  argument,  we  find  that  the  mean  cross-roll  wind  is  given 
approximately  by  V  =  2.5  ms~l .  The  distance  between  the  rolls  is  approximately  1 400  m. 
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which  produces  a  period  of  9.33  minutes  that  is  close  to  the  value  of  9.45  minutes  given 
by  (3.1).  We  next  check  these  values  for  the  period  against  the  nonlinear  results  for  two 
of  the  time-dependent  amplitude  coefficients  given  by  the  numerical  integration.  The 
period  given  by  each  graph  in  Figs.  3. 4-3. 5  is  approximately  9.3  minutes,  therefore 
demonstrating  that  the  model  integration  is  accurately  reproducing  the  expected 
periodicity  of  the  solution. 

3.2.  Case  Study 

We  use  the  physical  parameter  values  determined  in  Section  3. 1  to  examine  a  single 
case  of  forced  boundary  layer  roll  circulations.  Our  goal  is  to  determine  whether  the 
model  correctly  captures  the  patterns  produced  by  the  rolls  and  to  verify  the  behavior  of 
the  various  physical  quantities.  Because  we  have  chosen  a  case  in  which  the  solution  is 
damping,  the  magnitudes  of  the  physical  quantities  are  not  necessarily  representative  of  the 
actual  physical  solution.  However,  the  planview  patterns  and  profiles  produced  can  be 
analyzed  using  the  dimensionless  quantities  given  by  the  model. 

Our  first  step  is  to  determine  whether  the  model  captures  the  proper  spatial 
organization  of  the  roll  circulations.  One  way  of  examining  this  structure  is  to  analyze  the 
vertical  velocity  field  at  the  lower  boundary.  In  Fig.  3.6  we  present  such  an  analysis  in  the 
dimensionless  x*-y*~  plane  at  z*  -  0.  This  corresponds  to  the  height  of  the  lower 

boundary,  which  is  10  m  above  the  sea  surface.  The  areas  of  light  shading  indicate 
positive  vertical  velocities  and  the  areas  of  dark  shading  indicate  negative  vertical 
velocities.  There  are  distinct  bands  of  positive  and  negative  vertical  velocities  that  suggest 
the  model  is  capturing  the  spatial  organization  of  the  roll  circulations.  As  expected,  there 
is  no  net  vertical  motion  across  any  horizontal  plane. 


Fig.  3.6.  Planview  of  vertical  velocity  w'  in  the  dimensionless  coordinate  system  at  the  height  hiJs  of  the 
lower  boundary.  Areas  of  light  shading  indicate  positive  vertical  velocity  and  areas  of  dark  shading 
indicate  negative  vertical  velocity.  The  cellular  pattern  parallel  to  the  roll  is  an  artifact  of  the  graphics 
routine. 
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As  a  consequence  of  the  aspect  ratios  we  have  chosen,  the  roll  pattern  actually 
aligns  itself  a  bit  differently  in  the  dimensional  coordinate  system.  The  aspect  ratio  in  the 
x-direction,  ax  =  0.001,  corresponds  to  a  roll  wavelength  component  of  990  km  and  the 
aspect  ratio  in  they-direction,  ay  =  0.7,  corresponds  to  a  roll  wavelength  component  of 

1.414  km.  This  observation  means  that  in  dimensional  space,  the  rolls  would  actually  align 
themselves  almost  parallel  to  the  x-axis,  with  spacing  between  the  rolls  on  the  order  of 
1.4  km .  This  two-dimensional  planform  is  consistent  with  observations  of  roll  patterns, 
although  the  orientation  of  45°  to  the  right  of  the  mean  wind  is  not.  For  this  case  study, 
we  choose  such  a  large  orientation  angle  to  yield  a  periodic  solution  dominated  by  the 
advective  period. 

Another  way  of  examining  the  roll  structure  and  of  ensuring  that  the  rolls  are 
actually  transferring  heat  from  the  sea  surface  into  the  boundary  layer  is  to  analyze  the 
spatial  pattern  of  the  heat  flux  at  the  lower  boundary.  We  would  expect  areas  of  positive 
heat  flux  to  be  correlated  with  the  areas  of  positive  vertical  velocity.  In  Fig.  3.7  we 
present  this  vertical  heat  flux  analysis.  The  areas  of  light  shading  indicate  positive  heat 
flux  and  the  areas  of  dark  shading  indicate  negative  heat  flux.  As  can  be  seen,  the  areas  of 
positive  and  negative  heat  flux  are  well  correlated  with  the  areas  of  positive  and  negative 
vertical  velocity  in  Fig.  3.6. 

We  next  analyze  the  vertical  heat  flux  w'T  profile  integrated  over  the  x*~y* -plane. 
The  profile,  given  in  Fig.  3.8,  indicates  a  net  positive  heat  flux  over  the  entire  domain  that 
is  consistent  with  a  net  transfer  of  heat  from  the  sea  surface  into  the  boundary  layer.  We 
would  expect  this  heat  flux  to  be  largest  close  to  the  sea  surface  and  to  decrease  linearly 
with  height  as  indicated  by  observations  (Lilly  1966;  Stull  1988).  The  profile  in  Fig.  3.8  is 
consistent  with  this  linear  decrease  with  height  and  only  changes  slope  to  accommodate 
the  upper  boundary  condition  we  imposed  on  the  system.  The  modification  of 
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temperature  profile  by  the  roll  circulations  is  given  in  Fig.  3.9.  We  construct  this  profile 
using  temperature  amplitude  coefficients  and  Tm  only,  since  the  other 

contributions  would  integrate  to  zero  on  any  plane.  The  profile  suggests  that  the  strong 
positive  heat  flux  at  the  lower  boundary,  dominated  by  wavenumber  1,  does  indeed  change 
the  background  temperature  profile  T„{z)  the  most  there,  and  that  this  change  decreases 

linearly  with  height. 

The  fact  that  the  circulation  produces  a  linear  profile  means  that  the  observed  A ZT 
values  alone  cannot  be  used  to  estimate  the  correct  value  of  the  Rayleigh  number  Ra  for  a 
model.  Ideally,  we  would  set  Ra  =  0  in  the  model  and  compute  using  (2.29)  the  implied 
Ra  given  by  A  ZT  in  Fig.  3.9.  In  our  case,  the  effective  Ra  is  the  sum  of  the  one  imposed 
by  us  (2.75  x  105)  and  the  one  implied  by  Fig.  3.9  (4.50  x  103).  Here,  the  value  of  Ra  based 
on  observed  profiles  would  incorporate  the  effects  of  both,  with  each  one  reinforcing  the 
other. 

The  two  components  of  the  vertical  momentum  flux  profiles,  given  in  Figs.  3.10- 

з.  1 1,  are  found  by  integrating  uw*  or  vV*  over  the  x*-y*  -plane.  The  profile  given  by 
the  product  uw* ,  in  Fig.  3.10,  indicates  a  maximum  in  the  transfer  of  negative 
momentum  at  a  height  of  z  -  0.25,  which  corresponds  to  250  m  in  the  dimensional 
system.  This  suggests  that,  over  the  domain,  there  is  a  net  downward  transfer  of  westerly 
wind  u  to  the  lower  boundary  from  within  the  boundary  layer.  The  resulting  profile  for 

и,  given  in  Fig.  3. 12,  shows  that,  as  suggested,  momentum  has  been  transferred 
downward  giving  larger  negative  u  values  at  the  top  of  the  domain.  This  would  indicate 
that  the  circulations  should  tilt  with  height. 

The  profile  of  the  product  vV*,  given  in  Fig.  3.11,  indicates  a  maximum  in  the 
transfer  of  positive  momentum  at  150  m.  This  suggests  that  over  the  domain  there  is  a  net 
upward  transfer  of  southerly  wind  v*  from  the  lower  boundary  into  the  boundary  layer. 
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The  resulting  profile  for  v*,  given  in  Fig.  3.13,  shows  that  positive  momentum  is 
transferred  upward  giving  larger  positive  v*  values  at  the  top  of  the  domain.  The  shapes 
of  the  vertical  momentum  flux  profiles  and  associated  velocity  component  profiles  are  in 
agreement  with  the  results  found  by  B rummer  and  Busack  (1990). 

The  momentum  transfers  indicated  by  our  results  imply  that  the  initial  constant  wind 
profile  is  altered  to  give  a  southwest  wind  at  the  bottom  of  the  domain  and  a  more 
southerly  wind  at  the  top  of  the  domain.  Observations  would  reflect  the  sum  of  the 
original  wind  and  the  modified  wind,  meaning  that  an  observed  wind  turning  with  height 
might  originate  from  rolls  formed  in  an  environment  with  no  wind  shear.  These  effects, 
together  with  the  fact  that  wind  modification  is  predominantly  linear,  must  be  accounted 
for  when  using  observed  wind  profiles  as  model  input  (Haack  and  Shirer  1992). 

We  next  examine  the  vertical  velocity  profile  for  the  center  of  an  updraft  given  in 
Fig.  3. 14.  This  graph  indicates  that  in  the  updraft,  the  maximum  vertical  velocity  occurs  at 
the  lower  boundary  and  decreases  to  zero  at  the  top  of  the  domain.  This  result  is  not 
consistent  with  what  is  observed  in  the  atmospheric  boundary  layer.  We  would  expect 
from  observations  and  other  model  simulations  (e  g.  Schmidt  and  Schumann  1989)  that 
the  maximum  vertical  velocity  for  an  updraft  would  occur  in  the  middle  of  the  boundary 
layer  and  not  at  the  lower  boundary  (see  also  Figs.  2. 1-2.2). 

Therefore,  we  cannot  analyze  the  sea-surface  stress  patterns  given  by  mV*  and 
vV*  since  the  vertical  velocity  profiles  are  not  correct.  We  must  determine  why  this 
inconsistency  exists  and  correct  it  if  we  are  to  derive  useful  information  about  sea  surface 
stress  from  this  model.  In  the  next  chapter  we  discuss  some  possible  solutions  to  explore. 


/ 


(DIMENSIONLESS) 


Chapter  4 
CONCLUSIONS 


46 


The  case  study  results  in  this  thesis  show  that  the  model  does  in  fact  properly 
capture  the  spatial  organization  of  the  roll  circulations  and  does  correctly  reproduce  the 
vertical  profiles  of  heat  and  momentum  flux.  We  also  showed  that  the  temperature 
perturbation  is  a  maximum  at  the  lower  boundary,  just  as  we  would  expect  when  the 
primary  heating  source  is  located  at  the  lower  boundary.  This  result  is  in  good  agreement 
with  the  findings  of  Schmidt  and  Schumann  (1989)  that  are  based  on  a  study  of  a  higher 
resolution  LES  model.  Their  Fig.  18  illustrates  a  temperature  profile  that  does  indeed 
have  a  temperature  maximum  at  the  lower  boundary.  Our  results  also  reveal  that  the 
temporal  periodicity  of  the  system  is  consistent  with  both  physical  and  numerical  analyses 
of  the  solution  based  on  advection  by  the  cross-roll  mean  wind.  These  results  indicate  that 
the  system  of  equations  used  and  many  of  the  physical  assumptions  made  are  appropriate 
for  modeling  boundary  layer  rolls  in  a  marine  environment. 

Although  the  model  did  capture  many  of  the  expected  aspects  of  the  circulation, 
there  are  two  distinct  results  that  point  to  the  need  for  further  analysis  and  model  redesign. 
One  of  the  significant  results  of  the  model  integrations  is  that  we  could  not  find  a 
bifurcating  solution  that  came  into  a  state  of  equilibration.  As  mentioned  earlier,  there 
seems  to  be  a  source  of  energy  that  is  not  being  damped  properly  by  the  model.  Early 
investigation  points  to  excessive  growth  of  the  nonlinear  term  in  some  of  the 
thermodynamic  equations.  This  growth  manifests  itself  as  a  rapid  increase  in  the 
magnitude  of  the  temperature  amplitude  coefficient  Tm  that  represents  the  modification  of 
the  linear  temperature  profile  by  the  roll  circulations.  Although  in  the  case  of  the  damping 
solution  used  for  our  case  study  T0l  did  equilibrate,  it  was  to  a  value  that  greatly  exceeded 


47 


the  magnitudes  of  the  other  temperature  amplitude  coefficients.  This  mysterious  source  of 
energy  may  be  a  result  of  a  computer  coding  error,  although  we  performed  many 
integrations  using  symmetry  tests  to  identify  and  eliminate  this  type  of  error.  There  is  also 
a  chance  that  a  numerical  instability  causes  the  solution  to  grow  infinitely.  We  did  add  a 
numerical  dissipation  term,  as  described  in  Chapter  2,  that  seemed  to  control  the  growth 
and  so  allowed  us  to  at  least  find  a  damping  solution. 

The  other  major  concern  suggested  by  the  model  integrations  is  the  inconsistency 
in  the  velocity  field  at  the  lower  boundary.  The  vertical  velocity  profile  indicates  that  a 
maximum  in  upward  or  downward  motion  occurs  at  the  lower  boundary.  This  maximum 
however  should  occur  somewhere  near  the  midpoint  of  the  vertical  domain,  as  confirmed 
for  example  by  Fig.  18  of  Schmidt  and  Schubert  (1989)  and  Figs.  2. 1-2.2. 

There  are  two  ways  of  accepting  as  physically  plausible  the  results  of  the  present 
model,  although  neither  conform  to  the  original  concept.  First,  if  we  were  to  consider  the 
lower  boundary  to  be  near  the  midpoint  of  the  boundary  layer,  then  the  maximum  vertical 
velocity  occurring  below  it  would  be  consistent  with  observations.  Second,  if  we  limited 
our  domain  height  severely  to  be  within  the  lowest  part  of  the  boundary  layer  and  were 
interested  in  the  growth  of  the  circulation,  then  at  some  point  the  vertical  velocity 
maximum  would  occur  below  10  m  as  given  by  our  model  results. 

To  meet  our  original  intent,  however,  we  must  re-evaluate  our  application  of  the 
new  boundary  conditions,  which  not  only  use  nonzero  boundary  values  at  the  lower 
boundary,  but  rely  on  surface  layer  similarity  theory  to  specify  the  constants  in  them.  For 
example,  there  may  be  a  problem  in  how  the  values  of  the  lower  boundary  forcing 
constants,  $m  and  sT,  are  determined.  An  energetics  analysis  may  point  to  a  problem  in 
how  these  constants  contribute  to  the  boundary  source  and  sink  terms  that  result  from 
integration  by  parts.  The  more  likely  problem  is,  however,  that  the  way  in  which  we  apply 
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the  similarity  theory  does  not  result  in  a  proper  representation  of  vertical  fluxes  at  the 
lower  boundary. 

Upon  further  reflection,  it  seems  likely  that  the  similarity  theory,  which  is  used  to 
describe  the  dynamics  within  the  lowest  100  m  of  the  boundary  layer,  may  not  properly 
represent  the  velocity  scaling  we  need  for  specifying  the  lower  boundary  conditions  for 
kilometer-scale  circulations.  The  solutions  we  are  modeling  fill  the  entire  boundary  layer 
and  so  scale  with  boundary  layer  depth.  It  thus  seems  more  appropriate  to  use  convective 
scaling  to  specify  the  lower  boundary  conditions  for  velocity. 

Mathematically,  similarity  theory  gives  us  positive  values  for  the  two  constant 
forcing  parameters  when  the  boundary  layer  is  unstable.  This  forces  the  vertical  profiles 
for  temperature  and  vertical  velocity  to  take  the  form  of  the  vertical  basis  functions  given 
in  Figs.  3  . 1-3.2.  As  we  stated  earlier,  the  vertical  profile  for  temperature  is  consistent  with 
what  we  expect  at  the  lower  boundary.  A  maximum  of  temperature  at  the  lower  boundary 
follows  from  the  fact  that  T and  dT  jdz*  at  z  =  0  have  opposite  signs  when  sT  >  0.  This 
fact  can  be  seen  by  assuming  T ~  cos (z/zT)  at  z  =  10  m,  which  is  a  form  having  largest 

values  at  z*  =  0.  The  vertical  profile  for  vertical  velocity,  however,  does  not  allow  us  to 
obtain  a  maximum  above  the  height  of  the  lower  boundary  because  w*  and  dw*  / dz*  have 
opposite  signs  when  sm  >  0.  The  only  way  to  move  the  maximum  in  w*  into  the  domain  is 
to  reverse  the  sign  of  the  velocity  constant  forcing  parameter  sm  in  (2.35)-(2.37),  which  is 
easily  seen  by  assuming  w  ~  sin (z/zT)  at  z  =  10m  as  suggested  by  Figs.  2. 1-2.2. 

Hopefully,  we  can  solve  the  above  boundary  condition  problem  by  employing 
convective  scaling  theory  to  specify  the  lower  boundary  conditions  for  velocity  (e  g.  Stull 
1988).  We  hypothesize  that  this  scaling  theory  would  give  a  velocity  constant  forcing 
parameter  sm  that  is  negative  and  would  produce  vertical  basis  functions  that  yield  a 
maximum  in  vertical  velocity  above  the  lower  boundary.  Although  our  model  did  provide 
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an  accurate  representation  of  the  heat  and  momentum  flux  profiles,  we  must  recognize 
that  the  structure  of  these  profiles  may  be  correct  even  with  qualitatively  incorrect 
representations  of  vertical  velocity.  Until  the  boundary  conditions  are  correctly  specified 
to  produce  an  accurate  vertical  velocity  profile,  we  cannot  study  the  stress  variability 
pattern  at  the  sea  surface. 
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DETERMINING  THE  CONSTANTS  IN  THE  BOUNDARY  CONDITIONS 
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The  boundary  conditions  used  in  the  nonlinear  dynamic  model  are  based  on  the 
surface  layer  similarity  theory  from  Stull  (1988).  Two  empirical  convective  surface  layer 
relationships  are  used, 


kz  du 
«,  dz 


kz_dT^ 
T,  dz 


=  0.74(1 -94)  A=9.\t\. 


V 


(A.1) 

(A.  2) 


where  k  is  the  von  Karman  constant,  z  is  the  height  above  the  surface,  u  is  the  mean  roll- 
scale  wind  speed,  n.  is  the  friction  velocity,  L  is  the  Monin-Obukhov  length,  T  is  the 
mean  roll-scale  temperature,  and  T.  is  the  surface  layer  temperature  scale. 

The  boundary  conditions  were  developed  by  Ms.  Julie  L.  Schramm.  Her  goal  was 
to  make  the  boundary  conditions  (2.34)-(2.39)  consistent  with  the  similarity  relationships 
given  in  (A.  1)-(A.2).  The  relationship  in  (A.  1)  is  solved  for  u  giving  the  equation 


to--*. 

z. 


(A.  3) 
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where 


(A.  4) 


We  then  use  (A.  4)  and  its  derivative  to  eliminate  u,/k  to  give 


,  ]=du  [in (*/*.)-¥„(*/£)] 
dz{\ /z-^[¥m(z/L)]/^z} 


(A.5) 


Here  %  (z/L)  and  d[mm  (z/L)]/dz  are  obtained  from  (A.4)  using  the  definition  (A.  1)  of 
«9m  (z/L),  and  the  integral  evaluation  routine  from  the  DERIVE  software  package.  The 
equation  (A.5)  can  thus  be  expressed  in  the  form 


(A.6) 


We  next  introduce  the  relationships 


UK 

u  = - 

ZT  ~^LB 


(A.  7) 
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h LB  )  ^ LB  > 


(A.8) 


where  zT  is  the  height  of  the  top  of  the  boundary  layer  and  is  the  height  of  the  lower 
boundary.  Equation  (A. 6)  now  takes  the  form 

“V)  =  ^g-fm(z\zT^hLB,L)  (A.9) 

The  fimction  fm{z*,zT,z^hLB,L)  is  evaluated  at  z*  =  0,  which  corresponds  to  the  height 
z  =  hLB.  The  parameter  values  z„  and  L  are  determined  from  a  subroutine  supplied  by  Dr. 
George  S.  Young.  This  subroutine  was  implemented  by  Mr.  Dave  V.  Ledvina  during  his 
thesis  research  based  on  work  done  by  Dr.  Chris  W.  Fairall.  The  subroutine  uses  input 
values  for  the  mean  wind  speed,  humidity  difference,  and  temperature  difference  at  the 
lower  boundary  to  calculate  z0  and  L  The  evaluated  fimction  fm{ 0. zTlz0,hLB,L)  will  be 
referred  to  as  the  Schramm  momentum  constant,  sm,  which  now  represents  the  forcing  of 
the  wind  by  surface  layer  eddies  at  the  lower  boundary.  Example  subroutine  input  values, 
calculated  values  for  zo  and  L ,  and  associated  Schramm  momentum  constants  are  given 
in  Chapter  3. 

The  new  velocity  boundary  conditions  for  me  lower  boundary  that  apply  to  both  the 
u  and  v  components  of  the  wind  are  then 
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but  the  one  for  the  w  component  is  given  through  use  of  the  continuity  equation  by 


<^*(0)  <*V(o) 

dz*  m  dz ‘2 


=  0. 


(A.  12) 


Upper  boundary  conditions  must  be  applied  at  z*  =  1,  which  corresponds  to  the 
height  z  =  zT.  We  assume  a  rigid  stress-free  boundary  so  that 


du{  1)  =  dv*{\) 
dz  dz 


=  w*(l)  =  0. 


(A.  13) 


Boundary  conditions  for  temperature  are  found  with  a  similar  method.  The 
temperature  profile  given  in  (A.2)  may  be  solved  for  T  giving  the  relationship 


where 


(A.  14) 


(A.  15) 
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We  then  use  (A.  15)  and  its  derivative  to  eliminate  T./k  to  give 


7Y  ,  dT  [in (z/z0)-*PM(z/Z)] 
^  {\/z-^„(z/L)]/dz} 


(A.  16) 


As  with  velocity,  the  equation  can  now  be  expressed  in  the  form 


T(  z)  =  %fT(z’z°'L) 


(A.17) 


This  equation  is  now  made  dimensionless  using  the  relationships 


T_  TvkT„ 

S{ZT  -his) 


(A- 18) 


and  (A.8),  where  v  is  the  eddy  viscosity,  k  is  the  eddy  thermometric  conductivity,  and  T„ 
is  the  sea  surface  temperature.  Equation  (A.  1 7)  now  takes  the  form 


dT 

dz 


*  fr{Z  >ZTyZ°’hLB>L). 


(A.  19) 


The  function  fT[z* ,zT,z„hLB,  t)  is  evaluated  at  z  -  0  and  is  defined  as  the  Schramm 


temperature  constant,  sr,  which  represents  the  forcing  of  temperature  at  the  lower 
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boundary.  The  new  boundary  condition  for  the  lower  boundary  is  then 

r(°,  +  ,^U  (A.  20) 

We  assume  a  perfectly  conducting  upper  boundary  so  that 

r(l)  =  0.  (A.21) 

These  boundary  conditions  are  applied  to  the  Stokes  and  Helmholtz  problems  to 
determine  the  vertical  basis  functions  and  their  corresponding  eigenvalues,  as  summarized 
in  Appendix  B. 
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Appendix  B 

VERTICAL  BASIS  FUNCTIONS  AND  EIGENVALUES 

The  vertical  basis  function  and  corresponding  eigenvalues  for  temperature  are 
determined  by  solving  the  Helmholtz  problem 

rr  =  -V2r  (B.l) 

with  boundary  conditions  (A.20)-(A.21),  where 

<B2) 

T  is  a  dimensionless  temperature  eigenfunction,  and  y  is  the  corresponding  eigenvalue. 
We  use  separation  of  variables  to  find  solutions  of  the  form 


r  =trigl(x\y*)Fj(z'),  (B.3) 

where 

fr«o(*V)=l  (B4) 

trigx  ( x * ,  y* )  =  sin  2^fcc*  sin  2  nmy  (B .  5) 

trig2(x*,y *)  =  sin  2idx  coslnmy*  (B.6) 

trig3 (x\y')  -  cos2 nbc  sin  Irony  (B .  7) 


trigA(x‘ ,y')  =  cos2;rfx*  coslnmy*. 


(B.8) 


with  horizontal  wavenumbers  /  and  m,  and  F}(z *)  is  a  solution  of  the  eigenvalue  problem 

dz>i  +(r  a2)pj(z  )=° 

(B9) 

(B  10) 

^0)  =  o. 

(B  11) 

where 

o2  =  An2 a  2 12  +  4n2a2m2 
y-a2  =  Xj 

and  Aj  is  the  eigenvalue  of  the  vertical  basis  function  F}  (z*) .  We  number  these 

eigenvalues  by  increasing  value,  starting  with  A,  as  the  smallest. 

When  Aj>  0,  the  solution  is  of  the  form 

Fj  (z*)  =  A  sin 6>jz'  +Bcoso>jz*,  (B.  14) 

where  at2  =  kj ,  and  A  and  B  are  constant  coefficients.  Applying  the  first  boundary 


(B  12) 
(B  13) 
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condition  (B.  10)  produces 


B  +  sTa)]A  =  0.  (B.15) 


Applying  the  second  boundary  condition  (B.  11)  gives 

^sincj;  +  Bcoso)j =0.  (B.16) 


Solving  this  system  of  equations  yields 


Fj  (/)  =  cos a)fz*  - 


1 

sT<oj 


sin 


(B.17) 


where 


SjOHj  =  Xanatj.  (B.18) 

When  sT>l,  (B.14)  represents  the  only  solutions  that  exist.  When  sT  <  1,  however, 
the  first  eigenvalue  A,  is  negative.  This  eigenvalue  and  its  associated  vertical  basis 
function  Fx(z')  are  found  by  solving  the  differential  equation 


dz2 


(B-19) 
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r 


which  has  a  solution  of  the  form 

i: 


FHz)=Af*  (B.20) 

where  cy,2  =  A,  and  A  and  B  are  constant  coefficients.  Applying  the  boundary  conditions 
and  solving  the  resulting  pair  of  equations  gives 


where 


Fx{z*\  =  cosha),z* - sinh<y,z\ 

sTa  i 


(B  .21) 


sr<y,  =  tanhft), .  (B.22) 

The  vertical  basis  functions  and  corresponding  eigenvalues  for  velocity  are 
determined  by  solving  the  Stokes  problem 

aPax\  =  -PVV  +oxV/  (B.23) 

V  •  v*  =  0.  (B.24) 


with  boundary  conditions  (A.  10)-(A.  12),  where  P  is  the  Prandtl  number,  p*  is  the 
dimensionless  pressure,  ax  is  the  aspect  ratio  in  the  x  direction,  v*  is  a  dimensionless 
velocity  vector  eigenfunction,  and  a  is  the  corresponding  eigenvalue.  Equation  (B.24) 
together  with  the  boundary  conditions  imply  that  the  velocity  vector  v*  may  be  expressed 
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in  terms  of  the  curl  of  a  vector  streamfimction  A*,  this  vector  streamfunction  is  defined  as 
the  sum  of  two  vector  streamfunctions,  y/'  and  rf ,  where 


V* (x\yj,f)  =  ^£dV'„(t*)trigp(x\y)fiq(zy  +  Oj  +  Ok  (B  26) 

P= 0  q- 1 

n'(x\y\z\t *)  -  Oi  +  2Z  Wp(x*’^*K(z*)j  +  0k  (B 21) 

P= 0  4=1 

in  the  dimensionless  coordinate  system.  The  coefficients  y/pq(t‘)  and  r]pq(t *)  are  the 

temporally  dependent  amplitudes  for  the  pth  horizontally  dependent  trigonometric 
function  and  the  q\h  vertical  basis  function.  The  definitions  for  the  horizontally  dependent 
trigonometric  functions  are  given  by  (B.4)-(B.8).  The  function  hq(z‘)  is  the  qth  vertical 

basis  function  for  which  we  are  solving. 

Next,  we  curl  both  sides  of  the  Stckes  equation  to  eliminate  the  pressure: 


aVxv'  =  --V2(Vxv*l 

ax 

We  now  substitute  the  definition  for  v\  which  yields 

V2[V  x (V  x  A*)]  +axctV  x(Vx  A*)  =  0. 


(B.28) 


(B.29) 


We  substitute  for  A*  using  equations  (B.26)-(B.27),  after  applying  separation  of 


63 


variables,  and  expand  the  result  to  give  the  relationship 

‘^  +  (axa-a2)hq(z’)  =  °>  (B.30) 

where 

a2  =  4  n2a2l2  +  47z2a2m2  (B.3 1) 

axa  -a 2  =8q  (B.32) 

and  8  q  is  the  eigenvalue  of  the  vertical  basis  function  /^(z*). 

The  vertical  basis  function  and  associated  eigenvalues  are  now  determined  by 
solving  the  differential  equation 


+  ajq2hq(z)  =  0. 


where  mq  ~8q>  0,  and  by  using  the  boundary  conditions 


gW  ,g!W 

dz  m  dz2 


=  0 


W=0 


(B.33) 


(B.34) 


(B  35) 
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derived  in  Appendix  A.  In  this  case  we  have  arbitrarily  chosen  the  boundary  conditions 
for  the  velocity  component  w ,  although  we  could  have  chosen  either  of  the  other  two 
velocity  component  boundary  conditions.  The  solution  to  (B  33)  is  of  the  form 

hq(z*)  =  Csinmqz*  +  Dcosmqz\  (B.36) 

where  C  and  D  are  constant  coefficients.  Applying  boundary  condition  (B.34)  yields 

mqC  +  smtuq2D  =  0,  (B.37) 

and  applying  boundary  condition  (B.35)  gives 

Csinmq  +  Dcosmq  =0.  (B.38) 

Solving  the  pair  of  equations  (B.37)-(B.38)  results  in  the  vertical  basis  function 

hq(z*)  =  cosb79z*  -  smmq  sin  mqz\  (B.39) 

where 

~  COt  XU q.  (B.40) 

In  contrast  to  the  temperature  case,  these  are  the  only  solutions  that  exist,  regardless  of 
the  value  of  sm .  Thus  the  solutions  always  take  the  form  (B.36). 


